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ABSTRACT 

The luminosity and structure of neutron star magnetospheres are crucial to our 
understanding of pulsar and plerion emission. A solution found using the force-free 
approximation would be an interesting standard with which any model with more 
physics could be compared. Prior quasi-analytic force-free solutions may not be sta- 
ble, while prior time-dependent magnetohydrodynamic models used unphysical model 
parameters. We use a time-dependent relativistic force-free electrodynamics code with 
no free parameters to find a unique stationary solution for the axisymmetric rotat- 
ing pulsar magnetosphere in a Minkowski space-time in the case of no surface cur- 
re nts on the star. The solu tion is similar to the force-free quasi-analytic solution 
of ICont opoulos et al.l l)1999j) and the numerical magnetohydrodynamic solution of 
iKomissarovTij^obsj) . The magnetosphere structure and the usefulness of the classi- 
cal y-point in the general dissipative regime are discussed. The pulsar luminosity is 
found to be L « 0.99 ± 0.01/i 2 f2^/c 3 for dipole moment fi and angular frequency f2*. 
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1 INTRODUCTION 

Neutron star magnetospheres are suspec ted to have large re- 
gions of space that are nearly force-fr ee dGoldreich fc Julianl 
ll969l:lRuderman fc Sutherlandlll975l) . Self-consistent quasi- 
analytic solutions exist that describe the force-free (and 
some magnetohydrodynamic (MH D)) parts of the environ- 
ment of such systems (see, e.g ., | Contopoulos et all Il999t 
iGoodwin et alJl2004l lGruzinovll2005ir Quasi-analytic meth- 
ods cannot test the stability of their solutions and there is 
little hope to study the general nonaligned rotator. 

Recently. iKomissarovl i2005tl used a time-dependent nu- 
merical code to directly integrate the force-free and MHD 
equations of motion. They considered their force-free solu- 
tion unphysical due to uncontrolled fast reconnection in the 
current sheet that led to closed field lines far beyond the 
light cylinder. Their MHD version had no such defect, but 
they were forced to introduce unphysical evolution equations 
and parameters (see their section 4.2) in order to avoid nu- 
merical errors and strong unphysical features. They found 
a stationary MHD solution that may be unique, but they 
were uncertain whether their solution depended upon the 
unphysical model parameters. 

We investigate axisymmetric neutron star magneto- 
spheres by integrating the force-free equations of motion 
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using a newly developed general rel ativistic for c e-free e lec- 
trodynamics (GRFFE) extension jMcKinnevI l2005dl) to 
a general relativistic magnetohydrodynami cs (GRMHD) 
code called HARM dGammie et alJ 120031) . The origi- 
nal GRMHD code has been successfully used to study 
GRM HD models of accretion flows, winds, and jets 
jGammie et alJl200St iGammie Shapiro, fc McKirmevll2004l ; 
iMcKinnev fc Gammid 120041: iMcKinnevI l2005allbl<J) . The 
GRFFE formulation is a simple extension of HARM, so all 
the code developed for HARM is immediately useful and 
has identical convergence properties to the GRMHD ver- 
sion. We have shown that our GRFFE formulation is at 
least as accurate as the GRFFE formulation bv IKomissarovl 
J200ll2004l) iMcKinnevl2005dl) . Parabolic spatial interpola- 
tion and fourth-order Runge-Kutta tem poral integra t ion ar e 
used to improve accuracy. Compared to IKomissarovl |2005), 
our force-free scheme is improved by significantly lowering 
the reconnection rate in current sheets by forcing numeri- 
cally induced velocities into current sheets to vanish. This 
elimin ates all the difficulties encountered by IKomissarovl 
J2005D . 

Section|5]discusses the solution of the neutron star mag- 
netosphere. Section [3] s ummarizes the re sults of the paper. 
The notation follows lMisner et alJ il973l) and the signature 
of the metric is — h H — h Tensor components are given in a 
coordinate basis. The components of the tensors of interest 
are given by F^ v for the Faraday tensor, F for the dual 
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of the Faraday, and T M " for the stress-energy tensor. The 
field angular frequency is Of = F tr /F r ^ — Fte/Fe$. The 

* it 

magnetic field can be written as B 1 = F . The poloidal 
magnetospheric structure is defined by the (^-component 
of the vector potential (A$ = \&). The current system 
is defined by the current density (J) and the polar en- 
closed current {B^ = F^ t )- The el ectromagnetic l u minos - 
ity is L = -2tt LdBTTr 2 si n g. See iGammie et alJ <2003l) : 
iMcKinnev fc Gammiel ll2004ft ; lMcKinnevl (l2005dl) for details. 



2 NEUTRON STAR MAGNETOSPHERE 

Here the implementation details and the solution to the 
neutron star magnetosphere are discussed. It is most use- 
ful to the commun ity if a similar model to that chosen by 
iKomissarovl (I2005T) is studied and compared. This will show 
that one can study pulsar magnetospheres in the force-free 
limit using time-dependent numerical models. Des p ite th e 
unphysical parameters introduced by IKomissarovl {2005), 
our solutions are similar. 

The only additional note is that we and they choose a 
solution with no surface currents on the stellar surface by 
our and their choice of relaxed boundary conditions. This 
corresponds to, e.g., a neutron star with a crust in shear 
equilibrium (see, e.g.. lRuderman et alf l998). Other bound- 
ary conditions may lead to other solutions. 

2.1 Initial and Boundary Conditions 

As in lKomissarovl j2005l) . the magnetic field is approximated 
as an aligned dipolar field and the space-time is approxi- 
mated as Minkowski. The poloidal field at t — (and for all 
time for B r ) is set to be the perfect dipole with 

A ^ = ^Rl S m 2 8 = ^sm 2 9, (1) 
2r r 

where the magnetic dipole moment is fi — B po \ c Ri/2 for a 
polar mag netic field st r ength -B po i e and stellar radius R*. 

As in IKomissarovl j2005l) . the initial conditions are ar- 
bitrarily chosen to have a velocity of v l = and B^ = 0, 
and the proceeding violent non-stationary evolution even- 
tually relaxes to a steady-state solution . Thus, any solutio n 
found must be a stable solution. As in Komissarov (2005), 
the model is evolved for 55 light cylinder times. 

A steady-state solution with no discontinuities or sur- 
face currents on the stellar surface is found by choosing 
boundary conditions determi ned by an analysis of the Grad - 
Shafranov equation (see, e.g. . lBogovalovll997l : lBeskinll997l) . 
One is required to specify 2 constraints and to fix the 
magnetic fie l d com ponent perpendicular to the star. As in 
IKomissarovl d2005h . = F t<t> = 0, Q F = ft*, and B r are 
fixed in time, where S7* is the stellar angular velocity. This 
assumes the particle accelera tion gap is negligible, wh ich is 
not generally true (see, e.g.. iMestel fc Shibatal ll994'l. The 
velocity obeys th e frozen-in conditions (see equation 46 in 
lMcKirmevll2005dl) in st eady-state and axisymmetry. 

IKomissarovl i2005T) set the radius of the star to be 
Ri, — O.IRl, where Rl = c/Qf = c/O* is the light cylin- 
der. For i?* = 10 km, this means they chose a spin pe- 
riod of r » 2.1ms or O* ta 0.0207c 3 /GM for a neutron 
star with M = 1.44M©. We choose a similar frequency 



of O* « 0.0216c 3 /GM such that the light cylinder is at 
Rl = 46.3GM/c 3 . 

In practice, B r and Vlp = O* are fixed, while B 6 and 
Btf, are parabolically extrapolated into the neutron star sur- 
face from the computational domain. For a stationary, ax- 
isymmetric force- free solution, one can show that the field 
geom etry completely d etermines the velocity (see equation 
47 in lMcKinnevir2005dl) . For arbitrary field components B 1 
this prescription for if alwa ys leads to a time -like velocity 
within the light "cylinders" llMcKinnevll2005dl) . 

The polar axis boundary condition is such that the per- 
pendicular fluxes vanish. The outer boundary condition is 
obtained by extrapolating the field into the boundary zones 
and setting the velocity as described above, although the 
outer boundary is chosen to be far away to avoid reflections 
back onto the solution. In the event that there is a flux from 
the outer boundary into the computational grid, that flux is 
set to zero. 

2.2 Coordinates 

The computation is performed on a set of uniform rectangu- 
lar coordinates described by the vector field x^ , where each 
uniform coordinate is arbitrarily mapped to, e.g., t, r, 9, <j> in 
spherical polar coordinates. Apart from the code's ability 
to avoid significant numerical reconnection, the ability to 
choose an arbitrary 6 grid is crucial to obtain an accurate 
solution around the equatorial current sheet. 
The radial coordinate is chosen to be 

r = Ro + e x(1 \ (2) 

where Ro is chosen to concentrate the grid zones toward the 
surface (as Ro is increased from to J?*). An inner radius 
of Rin = R*, an o uter radius of R ou t = 2315GM/c 2 (50 
light cylinders as in IKomissarovl l2005h . and Ro = 0.9037?* 
are used. 

The 9 coordinate is chosen to be 

9 = tyx {2) + i(l - h{r)) sin(27T3; (2) ), (3) 

where labels an arbitrary uniform coordinate and h(r) 
is used to concentrate grid zones toward or away from the 
equator. The wide dead zone in inner-radial regions and the 
current sheet in outer-radial regions are both resolved using 

K r ) - | - + — atanf — ) ) (Tlouter — dinner) +^inaer, (4) 

\2 n \ r J J 

where h inncl = 1.8, h outcI = 0.05/i inllcr , r = 5GM/c 2 , and 
r s = 15GM/c 2 . 

Different resolutions were chosen to test convergence. 
In terms of radial vs. 9 resolution, we used 64 x 64, 80 x 64, 
80 x 128, 80 x 256, 80 x 128, 160 x 128, 160 x 256, and finally 
for the fiduc ial model we used 480 x 256 to reach a similar 
resolution of lKomissarovl (120051) . who used 496 x 244 for their 
final resolution. We find that all quantities have converged 
to within 1 — 5% at our highest resolution, where the exact 
percentage depends on the quantity. 

2.3 Magnetosphere Solution 

First, the code is checked that the solution reaches a steady- 
state. The solution does reach a well-defined steady-state by 
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Figure 1. Upper two lines show the value of Qp/Q* as a function 
of 9 for r = 0.2Rl (solid line) and r = 0.7Rl (dashed line) for 
64 9 zones and 80 radial zones. Lower lines show Qp/Q+ — 0.2 for 
256 9 zones and 48 radial zones. Dire ctly comparable to lower 
panels in figure 6 of lKomissarovl <2005lV 

about t ~ 1000GM/c 3 , or after about 20 light cylinder light 
crossing times, out to about r ~ 1000GM/c 3 for all of the 
domain except within the equatorial current sheet. The part 
of the equatorial current sheet that connects to the closed 
zone undergoes slow, small- amplitude oscillations. The fre- 
quency of these oscillations is proportional to the reconnec- 
tion rate allowed by our reconnection model, while the am- 
plitude is inversely proportional to the reconnection rate. 
Based upon the timing of the oscillations, the fiducial model 
has a reconnection period of ~ 10ms near the light cylinder. 
In the limit of an unlimited reconnection rate, there are no 
oscillations and the magnetosphere has a larger closed zone 
that extends some what beyond the light cylinder, similar to 
IKomissarovl (I20051) . 

Second, the code is checked for proper integration. 
Since the neutron star has f2* = Qf over the entire sur- 
face, then in axisymmetry and for a stationary flow, Qf 
should be the same for all field lines unless force-free (or 
ideal MHD) is violated. Notice that strong super-corotation 
(sub-corotation) would artificially enhance (diminish) the 
spindown luminosity and decrease (increase) the size of 
the closed zone. Figure shows Qp/Q* as a function of 
9 for r = Q . 2Rl a nd r = 0.7Rl, the same locations as in 
IKomissarovl J2005I) for their figure 6 (lower 2 panels). The 
plot is for t = 2546.3GM/c 3 w 55R L /c for a 9 resolution 
of only 64 zones and 80 radial zones (upper lines) and of 
256 6 zones by 480 ra dial zones (lower 2 lines displaced by 
0.2) . IKomissarovl (|2p05)'s figure shows that even at a 9 res- 
olution of 244 that the peak-to-peak fractional variation is 
22% near the equator and 45% overall (they have problems 
at the poles). Our model with only a 9 resolution of 64 has 
a peak-to-peak variation of only < 10%, and there are no 
significant artifacts at the poles. Over the entire domain, 
the spin of the magnetosphere deviates at worst by 5%. At 
a resolution of 480 x 256, the variation of Qf/^I* over the 



Figure 2. Flux function A$ in a box with size 3-Rl X 3Rl (Rl is 
light cylinder) as in figure 3 (upper left panel) of IKomissarovl 
(2005). There are 80 contours from A^, = (polar axis) to 
^0 = fi/R* = 9.57/if2*/c (equator on star). A single thick solid 
contour has been added for the closed field line that touches the 
theoretical light cylinder. The vertical solid line is the theoretical 
light cylinder Rl- 

computational grid out to r ~ 1000GA//c 3 (large radii still 
has initial transients) is < 3%. 

Finally, the region where the stationary solution is ac- 
tually force-free is checked, which is tracked by our code 
jMcKinnevll2005dl) . Force-free is only violated in the very 
thin current sheet and only in the sheet for r > 2.8Rl- En- 
ergy is not strictly conserved when force-free is violated, but 
the volume occupied by this region is negligible and so en- 
ergy loss is negligible and does not affect the luminosity as 
a function of radius described below. 

Figure [5] shows the magnetic flux function {A^). The 
theoretical light cylinder at R — Rl very closely follows 
the true outer light cylinder except in the very thin current 
sheet, where the light "cylinder" extends slightly radially 
outward by ~ 5%Rl- This magnctosph eric structure i s sim- 
ilar to the force- free solution of IContopoulos et alJ <ll999l) 
(although they have some kinks at the light cylinder) and 
the MHD solution of lKomissarovl (I2005T) . 

Figure |3] shows the structure of the magnetosphere near 
the slightly diss ipative current she et. In our solutions and 
the solutions of Komissarov (2005), there are always some 
very thin, extended closed field loops at the equator. Us- 
ing our low-resistivity model that keeps the current sheet 
from reconnecting, the heights of these extended closed field 
loops are smaller for increasing resolution, which also more 
sharply defines what is qualitatively associated with the so- 
called y-point. For general dissipative applications, however, 
the location of the y-point is ill-defined a nd cannot be quan- 
titatively determined. Notice that while IKomissarovl j2005l) 
state that their closed zone reaches all the way to the light 
cylinder, it is clear from their plots and other statements 
they make that there is no well-defined closed-open tran- 
sition, which is consistent with our results. However, one 
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Figure 3. Flux function for zoomed-in region of Fig- 

ure |5j around light cylinder. There are 80 contours from = 
1.09/^f2*/c to A^ = 1.80/jf!*/c. The right vertical line represents 
the theoretical light cylinder. The left vertical line points to the 
field kink point at the equator. 



can identify at lea st three quantitative measures, similar to 
IKomissarovl (I2005I) . 

First, the first closed field line to develop a kink at the 
equator near the light cylinder has an angle of 



9 kink « 76° ± 3° 



(5) 



between the field line and equator at the ki nk. This is simi- 
lar to the solution found bv lGruzinovl (|2005) . As they found, 
we find that this result is independent of f2* since the so- 
lution is smooth through the surface. For the model being 
presented, the intersection of the kinked field line with the 
stellar surface is at #kmk,* ~ 70° between a horizontal line 
and the field line at the stellar surface. The kink intersects 
the equator at 



J-kink « 43.08 ± QMGM/c = 0.931 ± 0.001i? L , 



(6) 



slightly within the theoretical light cylinder. The kink posi- 
tion oscillates in time by ~ ±0.01i?L, where the prior error 
estimate is smaller since a mean position in time was iso- 
lated. At the kink, the vector potential has a value 



1.35 ±0.01 



(7) 



which is quite close to the value of $ opcll = 1.36/xf2*/c 
given by IContopoulos et, all (1 1999}) and by t he va lue of 
tyy P = 1.375±0.005pf2*/c given bv lKomissarovl f2005) . That 
the solution kinks before the theoretical light cylinder could 
be due to the lack of a perfect dipole due to the light cylinder 
being close to the stellar surface, the remaining dissipation 
in the model, or could be intrinsic to our stable model. A 
comparison of the numerical dissipative y -point and the t he- 
oretical dissipationless y-point (see, e.g.. lUzdensk v 2003) is 
left for future work. 

The first kinked fi eld line appears ar ound r w 0.9Rl 
in the fiducial model of IKomissarovl i2005l) . and this agrees 



with our results. Since the nature of the dissipation is quite 
different in each code, this suggests that our similarly chosen 
fl* and so the lack of a perfect dipole, rather than dissipa- 
tion, is the reason why the first kinked field line appears 
inside the theoretical light cylinder. Otherwise perhaps it is 
intrinsic to an y stable solut i on. 

Second, in lKomissaro 3 J20051) they define the ' 'y-point" 
value of the vector potential to be at the theoretical light 
cylinder at 9 — tt/2, which we find to be 

c 



= flight = 1.27 ± 0.01 '- 



(8) 



which is quite close to the value of 'I'separat 
given b vlGruzinovl (I2005I). similar to <]/ 0] 
given m IContopoulos! J2005h . similar to ^ 
( 2003)7and similar to 



(2005) 



cj given bv IKomissarovl 
of the y-point value by IKomissarovl i200. 
which we defined as flight- One should com- 



1.27MV C 

open ' — ■ 

= 1.27/uf2*/c 
1.265 ± 
Notice that the 
is 



given by Timokhin 
0.005OO*/. 
definition 
their ^ yp 

pare our flight to their ^ yp , which are the same measure- 
ment. However, while they get 1.37, we get 1.27 in units of 
c = fi* = fl = 1. 

Third, as in IKomissarovl (l2005Tl . the vector potential 
value continues to drop along the equator until r > 5Rl 
where their value just oscillates around f opcn ~ 1.265 ± 
0.005(MVc) 



A 



4> | open 



We find that such a measurement gives 
.fJ.S~l* 



= 1.226 ±0.005^-, (9) 

c 

which i s quite close to the value of ^opon = 1.23/^fi*/c 
given in ICimtopflulflsj 12QQ5J ) and slightly different than the 
value given bv | Komissarovl 12005). That the value given by 
IKomissarovl J2005) is larger means they have more open flux. 
Thus, one expects their luminosity to be larger, which is the 
case. 

Of significant interest is the luminosity coefficient of the 
dipolar approximation. We find that the luminosity is 



0.99 ±0.01 



/'" 



(10) 



This result for L agrees with the result bv lGruzinovl l|20 
despite the detailed differences in the location of the kink- 
point, y-point, or the amount of open flux vs. closed flux. 
Since energy is explicitly conserved in our numerical method, 
the energy flux through each radius follows this same for- 
mula for the entire region that has rea ched a steady-state 
and does not violate force- free, unlike in IKomissarovl j2005l) 
where energy is not explicitly conserved. Note that for a 
lower resolution model, e.g. 160 radial zones by 128 9 zones, 
that L ~ 1.10 ± 0.05 in units wi th /j, = = c = 1 , which 
is quite similar to the result bv IKomissarovl i2005l) . Given 
that our code with 64 9 zones generated less error in Qf 
than even their model with 244 9 zones, then likely their 
solution with 496 x 244 radial vs. 9 zones is simply not as 
converged as our solution with 480 x 256 zones. However, 
the agreement between our force-free model and their ideal 
MHD model is impressive given the many unphysical pa- 
rameters they had to introduce to keep the solution realistic 
and the code stable (see their section 4.2). 

Figure [1] shows the large-scale structure of the field at 
the final time of 55Rl / c out to 30Rl ■ There are 40 contours 
from A<f, = to A<f, — 1.41/ jO^/c. The field i s essentially 
monopolar. This is similar to IKomissarovl i2005t) . although a 
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Figure 4. Flux function out to 30Rl ■ There are 40 contours 
from to 1.41/iS7*/c. 
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Figure 5. Current function Bj, vs. flux function \P = A^. 

reconnection event and plasmoid motion disrupts their field 
atr~ 23R L - 

A contour plot of the pulsar current function B<f, shows 
that the quantity follows field lines and has a significant 
current increase across the separatrix. A plot of the current 
structure agrees with the clas s ical m odel of pulsar current 
closure and with IKomissarovl (j200RT ) . To gain quantitative 
insight, Figure |K| shows the pulsar current function vs. the 
pulsar flux f unctio n, which is a similar plot to figure 4 in 
IKomissarovl (|2005h . We find that 

* ma x = 1.07^ at B = ±1.O5^, (11) 
at r = I.IRl- This is similar to the distribution shown by 



Contoooulos ct (Il999f) and is also similar to figure 4 in 
Komissarovl (l2005h . 



3 CONCLUSIONS 

A stationary force-free solution is found for the neutron star 
magnetosphere with a dipolar surface field in Minkowski 
space-time with no stellar surface current. The luminosity 
follows the standard dipolar luminosity with a coefficient 
determined accurately to be k — 0.9 9 ± 0.01. 

The difficulties encountered by IKomissarovl (120051) in 
the force-free regime were avoided, and the unphysical pa- 
rameters they introduced in the MHD regime were avoided. 
Their MHD solution is similar to our force-free solution, 
which suggests that the solution we and they find is accu- 
rate, stable, and may be unique. 
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